#ifndef _TPE_TERM_
#define _TPE_TERM_

#include "../DF_Common/basis_state.h"
#include "../DF_Common/const_nucl.h"
#include "../DF_Common/coupling_coeff.h"
#include "../DF_Common/progress_bar.h"
#include "../DF_Local/wavefun_coeff.h"
#include <boost/serialization/vector.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>

using std::cout;
using std::endl;
using std::vector;

struct t014_struct {
  int t_0;
  int t_1;
  int t_4;
};
struct l_struct {
  int l12p;
  int l12;
  int l3p;
  int l3;
};
struct Hc_struct {
  int lbar;
  int lbar_12;
  int lbar_3;
  int i_1;
  int i_2;
  int i_3;
  int i_4;
  double H0;
  double H1;
  double H2;
  double H3;
  double H4;
  double val;
};
struct Ic_struct {
  int lamb_1;
  int lamb_2;
  int lamb_3;
  int lamb_4;
  int lamb_i;
  int lamb_j;
  int lamb_k;
  int lamb_l;
  int lamb_12;
  int lamb_3k;
  int lamb_4l;
  int q_1;
  int q_2;
  int q_3;
  int q_4;
  double val;
  // int s_0;
  int s_1;
  int s_2;
  int s_3;
  int s_4;
};

struct Jaco_nlpnl_struct {
  // [nlp_num][channel][i]
  int i; // initial state index
  int f; // final state index
};

class Tpe_term {
public:
  Tpe_term(LEC_CONST &lec_t, Wave_Fun &wave_t, CG_coeff &cg_t,
           SixJ_coeff &sixJ_t, NinJ_coeff &ninJ_t, Jacobi &jaco_t)
      : LECs(lec_t), wave(wave_t), cg(cg_t), sixJ(sixJ_t), ninJ(ninJ_t),
        jaco(jaco_t) {}

  const LEC_CONST &LECs;
  const Wave_Fun &wave;
  const CG_coeff &cg;
  const SixJ_coeff &sixJ;
  const NinJ_coeff &ninJ;
  const Jacobi &jaco; /**/
  // initial
  int lbar123_max; // max of lbar lbar12 lbar3 calculate in Hc3_init
  JaJbJ JaJbJ_t;
  vector<double> cg_IntZero;
  int cg_J2max;

  int Ic1_size, Ic1_lbar_max;
  int Ic3_size, Ic3_lbar_max;
  int Ic4_size, Ic4_lbar_max;

  Binom bino;
  int bino_lmax;
  vector<vector<double>> BioHat;

  vector<double> Hat_vec;

  vector<vector<vector<vector<double>>>>
      T1S_sixJ;                    // [T'12][T12][S'12][S12] for T2 = 1
  vector<vector<double>> T3S_sixJ; //[S'12][S12] for T2=3 T'12=1 T12=1
  vector<vector<double>> T1_69J;   //[T12p][T12];
  double T3_69J;
  vector<vector<vector<double>>> S12_9J; //[S12p][S12][s_0]

  vector<l_struct> l_state;
  vector<vector<vector<vector<int>>>> l_vec; //[l12p][l12][l3p][l3]
  vector<t014_struct> t014_state;
  vector<vector<vector<int>>> t014_vec; // [t_0][t_1][t_4] = index

  vector<vector<t014_struct>> t014_state2;       //[l][i]
  vector<vector<vector<vector<int>>>> t014_vec2; // [l][t_0][t_1][t_4] = index

  vector<vector<Ic_struct>> Ic1_index;     // [l_state][index][15]
  vector<vector<double>> Ic1_val;          // [l_state][index]
  vector<vector<vector<double>>> Ic1_val2; //[l_state][index][14]

  vector<vector<vector<vector<Hc_struct>>>>
      Hc1_state; //[l_state][index][t014_state]
  vector<vector<vector<vector<vector<double>>>>>
      Hc1_val; //[l_state][index][t_index][index][H]

  vector<vector<Ic_struct>> Ic3_index;     // [l_state][index][15]
  vector<vector<double>> Ic3_val;          // [l_state][index]
  vector<vector<vector<double>>> Ic3_val2; //[l_state][index][14]

  vector<vector<vector<vector<Hc_struct>>>>
      Hc3_state; //[l_state][index][t014_state]
  vector<vector<vector<vector<vector<double>>>>>
      Hc3_val; //[l_state][index][t_index][index][H]

  vector<vector<vector<Ic_struct>>> Ic4_index;     // [s_0][l_state][index][15]
  vector<vector<vector<double>>> Ic4_val;          // [s_0][l_state][index]
  vector<vector<vector<vector<double>>>> Ic4_val2; // [s_0][l_state][index][14]

  vector<vector<vector<vector<vector<Hc_struct>>>>>
      Hc4_state; //[s_0][l_state][index][t014_state]
  vector<vector<vector<vector<vector<vector<double>>>>>>
      Hc4_val; //[s_0][l_state][index][t_index][index][H]
  vector<vector<vector<int>>> vec_flag_tot;
  //  TPE
  double pf_tpe_c1, pf_tpe_c3, pf_tpe_c4; // pre-facotr in TPE term

  vector<vector<vector<double>>> vec; // vec[channel][jaco_index a][jaco_index
                                      // b]
  vector<vector<vector<double>>> vec1, vec3, vec4;

  void init(vector<vector<int>> Jaco_l123_info, vector<int> Jaco_E_info,
            int jaco_J12_2max, int jaco_j3_2max);

  void init_cgInt(int Jmax);
  double cg_int(int ja, int jb, int J);

  void BinoHat_init(int lmax);
  double BinoHat_cal(int lbar, int i1) {
    if (lbar > bino_lmax) {
      cout << " wrong at cal_BinoHat lmax " << endl;
      exit(0);
    }
    return BioHat[lbar][i1];
  }

  void Ic1_init();
  void Hc1_init();

  void Ic3_init();
  void Hc3_init();

  void Ic4_init();
  void Hc4_init();

  int Hc1_cal(const int &l_index, const int &t0, const int &t1, const int &t4,
              const Ic_struct &Ic1, vector<Hc_struct> &Hc1);

  int Hc3_cal(const int &l_index, const int &t0, const int &t1, const int &t4,
              const Ic_struct &Ic3, vector<Hc_struct> &Hc3);

  int Hc4_cal(const int &l_index, const int &t0, const int &t1, const int &t4,
              const int &s_0, const Ic_struct &Ic4, vector<Hc_struct> &Hc4);

  double C1_term(const State_Jacobi &final_state,
                 const State_Jacobi &init_state);
  double C1_term2(const State_Jacobi &final_state,
                  const State_Jacobi &init_state);
  double C1_term3(const State_Jacobi &final_state,
                  const State_Jacobi &init_state);
  double C3_term(const State_Jacobi &final_state,
                 const State_Jacobi &init_state);
  double C3_term2(const State_Jacobi &final_state,
                  const State_Jacobi &init_state);
  double C4_term(const State_Jacobi &final_state,
                 const State_Jacobi &init_state);
  double C4_term2(const State_Jacobi &final_state,
                  const State_Jacobi &init_state);

  void build();

  void print(std::string filename);
  void print();
  void print(int ch_num);

private:
  int Tri(int &j2a, int &j2b, int &j2c);
  int Tri(const int &j2a, const int &j2b, const int &j2c) {
    int A = j2a;
    int B = j2b;
    int C = j2c;
    return Tri(A, B, C);
  }
  int Tri(int j2a, int j2b, int j2c, int unuse) { return Tri(j2a, j2b, j2c); }
  int phase(int n) { return n % 2 ? -1 : 1; }
};

#endif
